Stable and unstable attractors in Boolean networks 
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Boolean networks at the critical point have been a matter of debate for many years as, e.g., scaling 
of number of attractor with system size. Recently it was found that this number scales superpoly- 
nomially with system size, contrary to a common earlier expectation of sublinear scaling. We here 
point to the fact that these results are obtained using deterministic parallel update, where a large 
fraction of attractors in fact are an artifact of the updating scheme. This limits the significance of 
these results for biological systems where noise is omnipresent. We here take a fresh look at attrac- 
tors in Boolean networks with the original motivation of simplified models for biological systems in 
mind. We test stability of attractors w.r.t. infinitesimal deviations from synchronous update and 
find that most attractors found under parallel update are artifacts arising from the synchronous 
clocking mode. The remaining fraction of attractors are stable against fluctuating response delays. 
For this subset of stable attractors we observe sublinear scaling of the number of attractors with 
system size. 

PACS numbers: 89.75.Hc, 05.10.-a, 05.50.+q, 05.45.Xt 



C 

O 

(N 
> 

O 



O 

-i— > 



c 

o 
o 



X 



Boolean networks at the critical point (sometimes also 
called Kauffman networks) have been discussed as simpli- 
fied models for gene regulation networks for many years 
01 111 We currently experience a resurgence of inter- 
est in these models, as structure and dynamics of the 
genetic network in a living cell become visible thanks to 
new powerful experimental techniques (DNA chips) Q. 
From the theorist's point of view, Boolean networks ex- 
hibit interesting statistical mechanics with a prominent 
order/disorder phase transition |5|- Earlier, the critical 
state has been postulated to have some relevance in the 
biological context as the scaling properties of numbers of 
attractors with network size appeared to resemble how 
number of cell types scale with amount of genetic infor- 
mation when comparing simple and complex organisms 
Until recently it was believed that the total number 
of attractors increased as y/~N 0] . This has been falsified 
by improved simulation techniques [|| and it was shown 
that the total number of attractors grows faster than any 
polynomial 0,0. 

Let us here step back for a moment and reconsider 
Kauffman networks in the context of their original moti- 
vation, as models for biological systems. While the use 
of models discrete in time is an established approach in 
many circumstances of biological modeling, such ideal- 
izations always have to be treated with special care. In 
the case of Kauffman networks, the system evolves by 
a synchronous update of all nodes at integer values of 
time. Such a clocking, however, can produce spurious 
synchrony. For instance, subsystems are kept phase syn- 
chronized even if they are not interacting at all. In or- 
der to circumvent computational artifacts it has been 
suggested to use a more natural updating schedule @- 
For example it has been shown that the complex spatio- 



temporal patterns observed under synchronous update 
often disa ppe ar when units are updated asynchronously 

In this paper we reconsider Boolean networks at criti- 
cality, while destroying spurious synchrony by equipping 
the nodes with weakly fluctuating response delays. This 
allows us to analyze the stability of dynamical attractors 
in the discrete network model. A deterministic Kauff- 
man network at an attractor is perturbed by a slight 
shift of update events forward or backward in time. If all 
such perturbations die out, i.e. the system returns to the 
identical attractor, we call this attractor "stable" . Other- 
wise ongoing temporal fluctuations accumulate and drive 
the system away from the attractor. These latter cases 
correspond to attractors that are an artifact caused by 
synchronous update of the deterministic system. When 
systematically applying this method to Kauffman net- 
works we obtain as main result that the number of sta- 
ble attractors grows sublinear ly with system size (see Fig. 
Ma)). 

Let us study a Kauffman network composed of TV bi- 
nary nodes where each node determines its state Xi by 
applying a Boolean function (a rule table) fi : {0, 1} 2 — > 
{0, 1} on inputs received from two other nodes a(i) and 
b(i), according to a randomly chosen quenched topology. 
To be definite, we exclude self-couplings. Starting from 
an arbitrary initial condition (xi(0), #2(0), . . . , a;jv(0)), 
states of the nodes are synchronously updated at inte- 
ger times t according to the Boolean function 



Xi(t + 1) = fi(x a (i){t),X b (i)(t)). 



(1) 



The network itself as defined by fi, a(i), and b(i) remains 
constant in time. Running the system from a randomly 
chosen initial condition, its finite discrete state space of 
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FIG. 1: Frequency and accessibility of stable and unstable 
cycles in Kauffman networks, (a) Average number of stable 
(*) and unstable (□) attractors as a function of the number of 
nodes in the network, (b) Fraction of initial conditions leading 
to a stable attractor (solid line) and the ratio between num- 
bers of stable and all attractors (dashed line). Data points in 
(a) are averages over R independent realizations of the Kauff- 
man network, R = 10 7 for N < 24, R = 10 6 for 24 < N < 31 
and R = 10 4 otherwise. For increasing computational speed 



the networks are subject to the decimation procedure j(Sj be- 
fore simulation. For the decimated network we fully enumer- 
ate the state space such that it is certain that all cycles are 
detected. Sizes of basins of attraction in (b) have been es- 
timated in 10 5 networks, testing 100 randomly chosen initial 
conditions in each original network (no decimation applied). 



2 N possible states ensures that, eventually, a state reap- 
pears that has been encountered before. From thereon, 
the deterministic system will indefinitely follow the at- 
tractor it reached (which is either a periodic limit cycle 
or a fixed point). Different initial conditions may lead 
to the same or to a different attractor. The total num- 
ber of attractors is a characteristic property of a net- 
work. The expected number of cycles in an ensemble of 
random Kauffman networks has been shown to increase 
superpolynomially with system size TV 0. 

Let us now define a criterion for stability of an at- 
tractor in the presence of deviations from deterministic 
parallel update. For this purpose, we replace the dis- 
crete update times by a continuous time variable t where 
nodes may update at any point in time. Our goal is 



to slightly desynchronizc the dynamics of the network by 
shifting the individual updates of nodes to slightly earlier 
or later time points. To avoid that this generates spuri- 
ous spikes during transitional phases (e.g. when several 
signals interact that used to be simultaneous, but now 
arrive at slightly different times), nodes have to be pre- 
vented from switching on a time scale s much shorter 
than the original integer update time step (i.e., s <C 1)- 
This is implemented by a low pass filter that averages out 
fluctuations on short characteristic times scales s, namely 
by averaging over the input signal according to 



Xi(t + 1) = e 



(28)- 



t+s 



fi(Xa{i)(t),Xb(i)(t))dt - 1/2 



(2) 

where O is the step function with O(h) = 1 for h > 
and O(h) = otherwise. Let us briefly check how this 
works. Imagine node a(i) switches on at time t, node 
b(i) switches off at time t + e, while fa is the function 
and. Without low pass filter (s — 0), node i switches 
on at time t + 1 and off again at time t + 1 + e. When 
the switching time scale s of nodes is sufficiently slowed 
down, s > e, the spurious spike is filtered out, i.e. node i 
remains constant. Note that in the limit of fast switching 
time scale s — > 0, Eq. (|2J) converges towards Eq. Q). 
In particular, all synchronous solutions of Eq. (Q, i.e. 
solutions with nodes switching precisely at integer values 
of t, are solutions of Eq. J3J) for arbitrary s < 1/2, as 
well. 

Starting from such a synchronous attractor of a net- 
work, let us now perturb it at some time T by temporarily 
retarding parts of the switching events. Thereby, a subset 
of nodes that would normally change state at time T is 
kept frozen in their present state during the time interval 
[T,T + e] with e < s. After t = T + e, we let the sys- 
tem evolve as usual according to Eq. ©. Note that the 
original and the perturbed solution differ only on time 
intervals [t,t + e[ for integer t. In general, the time lag 
may propagate, i.e. for each integer t >T some units flip 
at time t while others flip at time t + e in the perturbed 
solution. If, however, there is a later integer time t* > T 
such that either no flips occur at time t* or no flips hap- 
pen at time t* + e, the perturbation has been overcome 
and the system has regained synchrony. We call an at- 
tractor stable if for all possible perturbations of the above 
type (i.e. for all possible permutations of perturbed and 
non-perturbed nodes) the system regains synchrony and 
the original attractor is stabilized within a finite time in- 
terval. Otherwise the attractor is called unstable. In real 
world situations with continuous noise, such unstable at- 
tractors will accumulate phase shifts that eventually shift 
the system into some other, stable attractor. 

With this we here choose a particularly simple crite- 
rion for the stability of an attractor in a Boolean network. 
The system is on a stable attractor if after each small 
perturbation it reaches the attractor again where, as a 
minimal perturbation, a small deceleration or accelera- 
tion of a switching event is used. On unstable attractors, 
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such time lags do not relax. Thus, ongoing perturba- 
tions eventually lead to a change in time ordering of the 
switching events and the system reaches a different at- 
tractor. Note that this scenario is much better suited 
as a stability criterion than stochastically adding or re- 
moving switching events [TH ITfil | , which does not allow 
for the limit of infinitesimally weak perturbations. The 
low pass filter characteristics used here is further moti- 
vated by the dynamics of biochemical switches [TtI ] where 
molecule concentrations typically respond slowly leading 
to an overall low-pass filter characteristics of the switch. 
Low pass filter characteristics is a natural property of 
genes and is a simple means of stabilization which is of 
low cost and ubiquitous in nature. We here model this 
mechanism as an effective time average over the input 
signal that suppresses short pulses. Note that under syn- 
chronous update of the model networks, pulses cannot 
be shorter than unit time such that the filter can be ne- 
glected. 

Applying the robustness criterion to random Boolean 
networks at criticality, one observes that the average 
number of all attractors, stable and unstable ones, grows 
much faster than the average number of stable attractors 
alone (see Fig. ^a)). In large networks, almost all at- 
tractors are expected to be unstable. Interestingly, the 
probability to reach a stable attractor from a random ini- 
tial condition is much larger than the fraction of stable 
attractors, as shown in Fig.^b). Thus, unstable attrac- 
tors typically have significantly smaller basins of attrac- 
tion than stable attractors. The main result is that, with 




rescaled system size (N/5) a 



FIG. 2: System size scaling of the number of stable attractors. 
Plotted as a function of the rescaled number of nodes (N/5) a 
with a = 0.3, 0.5, 0.7 (left to right). Error bars in (b) indicate 
standard deviation divided by y/R with ensemble size R (cf. 
Fig. [TJ . 

system size N, the number of stable attractors grows sub- 
linearly as ~ N a with a w 0.5, as shown in Fig. [5] A least 
squares fit of the form c + bN a fits best ((x 2 ) = 0.00013) 
with the parameter values a — 0.448, c = 1.107, and 
b = 0.108 (with a correlation coefficient for this fit of 
r = 0.999742). Further let us analyze the number I of 
states contained in the attractors. While stable attrac- 
tors are shorter on average than unstable attractors, the 
distribution of attractor lengths is broader for stable than 
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FIG. 3: Statistics of attractor lengths for networks with N = 
10 (thin curves) and N = 30 (thick curves). The cumulative 
distribution is shown for stable attractors (solid lines) and 
unstable attractors (dashed lines). For N = 10 nodes the 
average length of stable attractors is (Z) na t = 2.59, of unstable 
attractors (Z)art = 3.56; for TV = 30 we find (i) aa t = 5.50 and 

(Z)art = 7.16. 

for unstable attractors, as shown in Fig. [3J The major- 
ity of long attractors with lengths far above average are 
stable. 

Can we understand by a simple picture how unstable 
attractors differ from stable ones? Most obviously, un- 
stable attractors occur when the network falls into two or 
more non- interacting clusters. When all updates in one 
of the clusters are delayed by the time e, this phase lag 
with respect to all other clusters cannot heal. All attrac- 
tors with flipping events in more than one network cluster 
are unstable. However, also in networks consisting of a 
single cluster (more precise: with a single strongly con- 
nected component) unstable attractors are found. Figure 
0] illustrates the coexistence of a stable and an unstable 
attractor in a small connected network. The example 
suggests that an attractor is stable if there is a single 
cascade of switching events. Let us consider the minimal 
number of simultaneous flipping events 

m = min\{i\x i (t)^x i (t+l)}\ (3) 

for a given attractor. The attractors with m = are 
the fixed points. These are stable by definition because 
no flipping events are to be retarded. Attractors with 
m = 1 are stable as well. These attractors contain a step 
with only a single flipping event. Going through this 
step the system always regains synchrony. For m > 2 
the attractor is likely to contain several chains of causal 
events as in Fig. 0fc) . In the simulations we find that a 
large fraction u of the attractors with m > 2 is unsta- 
ble, u = 0.856, 0.882, 0.899, 0.9094 for N = 10, 20, 30, 36, 
respectively. Thus the minimal number of simultaneous 
flipping events m allows for an almost perfect distinction 
between stable and unstable attractors. Note that m is 
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(a) 




(b) 
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FIG. 4: A stable and an unstable attractor in the same 
system, (a) Three nodes forming a directed cycle (feedback 
loop). Each node i performs the Boolean function negation on 
the input from its predecessor j, i.e. Xi(t+1) = 1 — Xj(t). (b) 
Temporal evolution of a stable attractor. There is a unique 
causal chain of flipping events (thick arrows) . A retarded up- 
date (shaded box) retards all subsequent events by the same 
amount of time. Thus perturbations heal immediately, (c) 
Unstable attractor. It can be interpreted as the superposi- 
tion of three independent chains of flipping events propagat- 
ing along the cycle of nodes. One of the chains is indicated by 
thick arrows as in (b). Retarding an event effects subsequent 
events in the same causal chain only. Causal chains remain 
phase shifted. The system does not regain synchrony after 
a perturbation. See reference fl8| for a detailed analysis of 
attractor stability in small systems. 



measured in the decimated networks |6J. 

Comparing these results to past studies of random 
Boolean networks at criticality (K ~ 2 inputs per node) , 
we obtain a distinctly different picture: Only a small 
fraction of all attractors are at all stable against small 
amounts of noise. Or, put differently, the effect of spu- 
rious synchronization due to a parallel update mode has 
been underestimated in previous studies, at least where 
these studies have been made with a potential applica- 
tion to biological systems in mind. In particular, charac- 
teristic properties of the attractor statistics are different 
when considering the subset of stable attractors: The av- 
erage number of stable attractors scales less than linearly 
while the number of unstable attractors shows a faster, 
superlinear growth with N. Also, stable attractors have 
a significantly larger basin of attraction than unstable 
ones. One may speculate that this latter property might 
have been the reason for the long prevalence of the opin- 
ion that the total number of attractors scales as y/N 0] • 
Mainly these stable attractors were likely to be found in 
the early studies using sparse sampling. 

If one aims at discussing Boolean networks as simple 
models for biological systems, our study suggests to con- 
sider more carefully the question of which attractors are 
at all relevant to the biological system. For example, 
Kauffman's observation of the number of attractors in 
critical random Boolean networks exhibiting a similar 
scaling with system size as the number of cell types with 
genome size in organisms Q seems to be wrong in the 
light of the results by Troein and Samuelsson |2J . How- 
ever, it appears to be still open to debate when consid- 
ering solely the subset of stable attractors. 
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